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On the dynamics of approximating 
schemes for dissipative nonlinear equations 


By Don A. Jones 


1. Motivation and objectives 

Since one can rarely write down the analytical solutions to nonlinear dissipative 
partial differential equations (PDEs), it is important to understand whether, and in 
what sense, the behavior of approximating schemes to these equations reflects the 
true dynamics of the original equations. Further, because standard error estimates 
between approximations of the true solutions coming from spectral methods — finite 
difference or finite element schemes, for example— and the exact solutions grow 
exponentially in time, this analysis provides little value in understanding the infinite 
time behavior of a given approximating scheme. 

The notion of the global attractor has been useful in quantifying the infinite time 
behavior of dissipative PDEs, such as the Navier-Stokes equations. Loosely speak- 
ing, the global attractor is all that remains of a sufficiently large bounded set in 
phase space mapped infinitely forward in time under the evolution of the PDE. 
Though the attractor has been shown to have some nice properties — it is compact, 
connected, and finite dimensional, for example— it is in general quite complicated 

Nevertheless, the global attractor gives a way to understand how the (infinite 
time) behavior of approximating schemes such as the ones coming from a finite 
difference, finite element, or spectral method relates to that of the original PDE. 
Indeed, one can often show that such approximations also have a global attractor' 
We therefore only need to understand how the structure of the attractor for the 
PDE behaves under approximation. This is by no means a trivial task. Several 
interesting results have been obtained in this direction. However, we will not go 
into the details. We mention here that approximations generally lose information 
about the system no matter how accurate they are. There are examples that show 
certain parts of the attractor may be lost by arbitrary small perturbations of the 
original equations (see Humphries, Jones, Stuart, 1993, and the references therein 
for a description of some of the results). 

Under certain hypothesis on the approximation, one can be guaranteed some 
structures of the attractor survive the approximation. For example, consider the 
ordinary differential equation (ODE) 

dx 

Tt =X{x) ' (i) 

where we suppose that x £ R n and that X is a C l function. Further, we suppose 
that the system (1) is dissipative and hence has a global attractor. Now suppose 
that (1) is approximated by 

§ = *<*) + >'(„), (2) 
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where y € JR" and Y is a C 1 function. Suppose further that 

imic* <* 

for some suitable e > 0. That is, (1) and (2) may be viewed as small C 1 per- 
turbations of one another. This seems to be a natural conditmn to reqmre of a 
perturbation in order to say something about how the global attractor of (2) rela 
to that of (1). Indeed, such systems have been studied by several authors and im 
creasingly stronger results have been obtained (see e.g. Phss & Sell, 1992, and 
references therefn). It is known, for example, that normally hyperbolic, invariant 

manifolds persist under such perturbations. 

to order to apply these results to PDEs, one must first construct firnte ^systems 
of ODEs that have the same global attractor as the infinite-dimensional PDE. T 
has been done for several dissipative PDEs including, for example, the Kuramoto- 
Sivashinsky equation, Cahn-Hilliard equation, Ginzburg-Landau certain reaction- 
diffusion equations, and the Navier-Stokes equations, Kwak, 1991. Such systems 

ao-e called inertial forms. f -i 

To be more specific, each of these PDEs can be viewed as an ordinary differential 

equation on a suitably chosen Hilbert space, H. We denote by ( , •) the tuner product 

and | - | the norm on H. Then these equations take the form 


du 

dt 


+ Au + R(u) = / 


( 3 ) 


u(0) = uo- 

Typically, the operator Au is -V 2 with Dirichlet or periodic boundary conditions. 
For the Navier-Stokes equations, for example, the term R{u) is the divergence ree 

part of (u • V)u (see Temam 1988). r TF ', 

In all cases but the Navier-Stokes equations, the existence of inertial forms (IF 
has been proven by showing the existence of an Inertial Manifold^ To date , inertid 
manifolds have been constructed as a graph m phase space of a Lipschitz function 
n“as, Sell, Temam, 1988). An inertial manifold (IM) for a dissipative evo- 
lution partial differential equation is a smooth finite-dimensional manifold ra p ase 
space, which is positively invariant under the solution operator and which umforady 
attracts every bounded subset of phase space at an exponential rate. It is clear 
if the IM exists, then it must contain the global attractor. Moreover, the reduction 
of the partial differential equation to the IM yields the inertial form. 

We denote by P the orthogonal projection of the space H onto the span of the 
first M eigenfunctions of A, and Q = I - P- We set p = Pu,q = Qu. Then the 
evolution equation (3) is equivalent to the system 
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If the IM is given as a graph of a Lipschitz function : PH i— ► QD(A) for 
M sufficiently large, then on this manifold the solutions of (3) are of the form 
u(t) = p(t) + <I>(p(tf)). Moreover, in this case, the inertial form is given by 

^ + Ap + PR(p + $(p)) =Pf P € PH. (4) 

Equation (4) does not have the same solutions as Equation (3) (only on M ). 
Rather, it has the same infinite time behavior as original PDE. Most importantly, 
it is an ODE. In view of the work of Pliss & Sell (1992) mentioned above, there is 
an advantage in approximating (4) in the C 1 sense. A candidate for such a system 
is 

^ + Ap + PR(p + $app(p)) = Pf p e PH, (5) 

with 


sup (|4($ app (p) - $(p))| + ||A(D$ aJ) p(p) - i?^(p))||£(i>//,QH)) < e, 
p€PH 


and where D'if denotes the Frechet derivative of the function \£. In this case, under 
reasonable assumptions on /?, the vector field in the approximate inertial form (5) 
may be viewed as a small C 1 perturbation of the vector field in the inertial form 

( 4 ). 


2. Accomplishments 

The main goal of this method of reduction is to implement the reduced ordinary 
differential system (5) in long-time simulations of solutions to the PDE, (3). Even 
in the case that the IM or a smooth function $ does not exist, the theory suggests 
looking for a global function $ a pp whose graph in phase space approximates the 
attractor. Indeed, many $app have been constructed. These approximations have 
been implemented in numerical schemes for a variety of equations and settings 
(see Jones, Margolin, Titi, 1993, and the references therein). We will discuss the 
effectiveness of these schemes below. 

Perhaps the most important role the IF, Equation (4), can play, as mentioned 
above, is to understand how the dynamics of approximating schemes relates to that 
of the original PDE. The first attempt at approximating $ in the C 1 sense was 
in Jones, Titi, 1993. There, $ was viewed as the asymptotically stable stationary 
solution of a certain PDE. One can then approximate $ by integrating this PDE 
forward for a short time. 

However, the situation may be much simpler than this. Consider the spectral 
approximation of (3) based on the eigenfunction of the linear operator A. One 
obtains the approximation 

du n 

—r— + Aa/v + PnR(um) = P^f 


( 6 ) 
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with initial data ii/v(0) = uo,jv* As shown in Foias, Sell, Temani, 1988, and Foias, 
Sell, Titi, 1988, if N is chosen sufficiently large, there exists a global function 
such that Mjv =Graph($/v) is an inertial manifold for (6). On this manifold, 
solutions are of the form u n ( t) = pN(t) + *& N{pN{t)) with pN(t) — Ptx/v(i). Further, 
on this manifold, (6) reduces to 

- + Apjv + PR(pn + $aKj>n(<))) = Pf, (?) 

at 

where P is defined above. Notice also that Equation (7) remains of dimension M, 
the same dimension as (4), a s N — ► oo. Again, it is not that (7) has the same 
solutions as (6) (only on Mn), but rather it has the same global attractor as (6) 
(since the IM contains the attractor). 

Now Equation (7) will play the role of the approximate inertial form, Equation 
(4). That is, we take $ apfl = tk/v- Moreover, it was shown in these two papers that 

c 

sup |A($(p) - <Mp)I < 

p€PH \ N 

for some 0 < 0 < 1/2 which depends on the nonlinear term R. However, more is 
true. 

Theorem Suppose that M is so large (determined by the spectral properties of A) 
that $, as described above exist. Then for all e > 0 there exists a N(e) such that 

sup [|A($(p) - $n(p))| + ||A(D$(p) - I , $n(p))||c(pw,c?<)] ^ e 

p€PH 

for all N > N(e). 

Proof. See Jones, Titi (1993). 

The IF, Equation (4), has the same dynamics as the original PDE. Moreover, 
the above shows that the spectral method based on the eigenfunctions of A is a 
small C 1 perturbation of the IF for N sufficiently large, since Equation (6) and (7) 
have the same attractor. Thus, this spectral method preserves certain structures of 
the attractor of the PDE, for example, the ones studied in Pliss, Sell, 1991, for N 
sufficiently large. 

A similar type of analysis may be possible for finite element methods. To do this 
properly, we should turn to a specific PDE. However, we will attempt to keep the 
exposition as general as possible. We denote by {V*}a>o u finite dimensional sub- 
space of differentiable functions (most typically piecewise linear functions), where 
one can think of ft as being the maximum partition size. Then one attempts to 
approximate solutions u(t) of (3) by functions u*(t) in V*. The functions u h solve 

(«? , X) + (A 1 ' V, A j / 2 X ) + (R(u h ), X ) = (/, X) 

u A ( 0 ) = 4 e v\ 


( 8 ) 
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where one can think of A 1 / 2 as d/dx in the ID case. 

The operator A 1 * can be defined from the equation (.4 A j/>, = ( .4 1 / 2 0 A l / 2 \) 

Further, projecting R and / onto the space V h , Equation (8) takes the form 

dxi ^ 

— + A h u h + R\u k ) = f k . (9) 

Since A is assumed to be self-adjoint, A h is also. Moreover, the spectrum of A h 
can be shown to approximate that of A. Thus, the space V h may be decomposed 

c n~ P - V 0 @ V ' In the same manner that the $ N was constructed in Foias, 
Sell, Titi, 1988, a global function $ A may be constructed for Equation (9) (see Jones’ 
Stuart, 1993, for the details) such that M h =Graph(4> A ) is an inertial manifold for 
(9). On this manifold, solutions are of the form = p h (t) + where 

p = P h u h . On the IM, (9) reduces to 

+ Ahph + phRh ( ph + = p k f do) 

for h sufficiently small. As in the case of the spectral method, the dynamics of (10) 
are the same as that of (9). Moreover, the dimension of (10) remains fixed (roughly 
on the order of M) as h -*■ 0. One would like to show C 1 closeness of (10) and (4). 
However, at this point all we have is the following 

Theorem For h sufficiently small there exists a function $ A such that (10) holds. 
Moreover, 

(i) for any p G PH there exists C(p) > 0 such that 

II (P + *(P)) - (P h p + $ A (P A p))|| < C{p)h ; 

(ii) for any p h G P h H there exists C(p h ) > 0 such that 

||(Pp A + $(Pp A )) - (p h + $ A (p A ))|| < C{p h )h. 


Proof. See Jones, Stuart (1993). 

On the Practical Side 

The above theory suggests that there may be an advantage in enslaving the high 
Fourier modes (in the case of the spectral method based on the eigenfunctions 
of A) in terms of the lower modes through the function $ app . Shortly after the 
discovery of the IF, Equation (4), many $ app were constructed and studied for 
various equations (see Jones, Margolin, Titi, 1993, and the references therein), 
including the Navier-Stokes equations (see for example Jolly, 1993). Schemes based 
on enslaving q ss <t>a PP (p) are generally referred to as nonlinear Galerkin methods 
since $ app = 0 gives the standard Galerkin scheme. 
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In Jones Margolin, Titi, 1993, we evaluate the effectiveness of the nonlinear 
Galerkin metho/in the context of spectral method (such schemes have now been 
constructed for finite element and finite difference schemes; see the references in 
Jones it et al. 1993). The goal of this work is to understand under what conditions 
the nonlinear Galerkin methods lead to a significant improvement in accuracy over 
the standard Galerkin method from a purely numerical analysis point of view. 

Recall that in general, if one approximates a smooth function u with respect 
to sor^e basis, the rate of convergence is limited by the smoothness of the basis 
functions. If the basis elements are C°° functions as in the case of the eigenfunction 
of the Unear operator A, the rate of convergence is only limited by the smoothne 
of u and compatibility of the function n with the basis elements of the expansion 
at the boundary (the presence of Gibb’s phenomenon, for example). 

Thus, if the solutions u(t) of equation (3) are very regular and compatible with 
the eigenfunctions {<,,} of A at the boundary, then the Fourier coefficients of the 

l ,i;fL mav decav very rapidly in wave number. Indeed, Foias, Temam, 1989, 
showed that, under such circumstances, solutions of the Navier-Stokes equation may 
decay exponentially in Fourier space. Similar results hold for the other equations 
mentioned above. Thus, the business of trying to approximate the q part of so utions 
via the function $ apP may not be effective when the q = Qu part oi f t he ^lu ‘ 
is exponentially small. That is, the approximation - 0, which leads to 

standard Galerkin scheme, may already be good enough. It turns out t ^ at W ^ 
controls the regularity and compatibility of solutions coming from the NSE and 
related equations is the compatibility with the basis functions and regularity of the 

f °We 6 dlmon{trate this for the Kuramoto-Sivashinsky equation (KS). A similar 
analysis holds for the NSE. This equation is given by 


dt dx 4 dx 2 


dx 


u( 0,x) = u 0 (r) 

u(t, x) = u(t, x + L) L > 0, t > 0. 

The KS equations appears in physics literature with f = 0. Here we b iave added an 
additional feature to the KS equation, namely, a forcing term /. We will use this 
forcing term to control the level of regularity of the solutions to the KS equation. 
It is clear that whenever /(*) is an odd function then the space of odd functions is 
invariant under the solution operator for the KS equation. For simplicity, we will 
restrict ourselves to the odd case. Hence, under these assumptions, one can easily 
show that the KS equation is equivalent to the evolution equation 


du 

dt 


+ Au — A l ^ 2 u + B(u,u) — /, 


u(0) = wo 
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on the Hilbert space H - {u £ Z 2 ((0, L))\u{x) = u(x + L),u(x) = ~u(L - x),x e 

St}. Here A — 3r « , the eigenvalues of A are A m = (2nm/L) i corresponding to the 
eigenfunctions <p m = sin(2nmx/L), for m = 1, 2, .. .. 

In either case^one can define the operators .4° for a > 0. One defines D(A a ) — 
{u € H\u = £lj=i u i‘r>hYl°jLi ^ 2 j a \uj\ 2 < oo}. Consequently, functions in D(A a ) 
are more regular, and more compatible (whenever it applies) with the eigenfunctions 
of A at the boundary, for larger a. 

We may then approximate solutions of the KS equation with either the standard 
Galerkin scheme 


+Ayn+A'' 2 y n + P n B(y n ,y n ) = P n f ( 11 ) 

or by the nonlinear Galerkin method 


dt + Az n + A 1/2 Z n + PnB(z n +$ app ( Zn ) >Zn + $ app ( 2n )) - p n f ( 12 ) 
for some clever choice for $ app . 

As mentioned, the rate of convergence of the two schemes is tied to the smoothness 
of the solutions and compatibility of the solutions with the basis functions (Since we 
are considering the KS equation with periodic boundary conditions, compatibility 
of the solutions is not an issue here. However, such cases are studied in Jones, 
Margolin, Titi, 1993 .) This is, in turn, tied to the smoothness and compatibility 
of the forcing term /. The rates of convergence of the two schemes is given by the 
following two theorems which is based on the work of Devulder, Marion, Titi, 1993 , 
and whose proofs can be found in Jones, Margolin, Titi, 1993 . 

Theorem Let u(t) = p(t ) + 9 (t) be a solution of the KS equation with u 0 on the 
attractor and f G D(A a ). Suppose y n solves (11) with y„(0) = P n u(0). Then 

M0-y»(0IU*<5r£r. 

A n+ 1 


In general, requiring / 6 D(A°) for some a > 0 requires not only that / be 
smooth, but also that / and its derivatives up to order 4a must satisfy the boundary 
conditions. Now we suppose that $ app satisfies certain conditions described in 

Jones, Margolin, Titi, 1993— such $ app abound. For the nonlinear Galerkin method 
we have ’ 

T heo r em Let u(t) be as in the previous theorem . Suppose z n solves (12) with 
z n { 0) = P n u( 0). Then [ ' 


I w (0 ( z n(t) + ^app(^n(0))IU 2 — Ta+1 


A n+ 1 
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for some 0 > 1 . 

For most $ aDD , 0 is not larger than two. Now one can see the issue here. As / 
becomes more smooth, larger a, the difference in the theoretical rates of convergence 

nf the two schemes decreases. , - , 

Let us push the smoothness of the solutions to an extreme. Consider the forced 

KS equation and suppose that for some a > 0, / € D(e ), t at is, 

Gevrey class (real analytic). Notice that 


>=1 


*'\‘ /4 |/>| 2 <o°, 


where f = T°° , /,¥>,. Under these assumptions for the forced KS equation, we 
lie (Propo^itTon 3.6 of Jones, Margolin, Titi, 1993) that the Fourier expansion 
of the solution converges exponentially fast. This means that the so utions are 
infinitely compatible with the basis functions at the boundary and analytic inside 
he domain. Hence, the high Fourier modes of the solutions have exponentially 
small norms, and there may be little advantage in approximating them. Indeed, in 

this case, we have 


_ r 1/4 

0 -< 7 l A m + 1 


IK*) - ym(<)lk J ^ Ci (*)“XhT 


where y m solves the Galerkin scheme (11) 
As in the above theorem, we have 


||u(t) - (Zm(f) + $app(Zm(<)))IU J ^ W)' 


_ vW 4 
B — ^1 A m + 1 

A m+ 1 


where z m is the solution of the nonlinear Galerkin scheme (12). 

Thus in the case that the solutions have Gevrey class regularity (spatially real 
analytic) the nonlinear Galerkin method only leads to algebraic improvements in 
the upper bounds of the rates of convergence over the standard Galerkin scheme. 
This little improvement might not be significant in computations. Of course, the 
overall improvement depends on how small the constant <7! is. 

We demonstrate this numerically. We choose $ apP to be 

(p) = A~ l Qm{f ~ S(p,p)), 

which we first studied for the NSE in Foias, Manley, Temam, 1988. For this choice, 
0 = 7/4. For our first example, we force the KS equation with 
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Notice that for this choice / 6 D(A a ) for a < 1/8. 

We have found a stable periodic orbit with this forcing. Such trajectories are on 
the attractor. To obtain the “exact” periodic solution, we run a Galerkin scheme 
using 100 modes; when we reach .5 time units, we start recording the data every 
.002 time units up to a time of .7 time units (Figure 1). If the initial data is taken 
near the periodic orbit, it will take some .4 time units to converge to the periodic 
solution. Since the periodic solution is on the attractor, we expect that the rate of 
convergence in this case to be of the form 

E gal = max ||u(f) - y m (t)\\it < ~~ , 

A m+ 1 

of course, here we are ignoring errors due to the discretization of time. Also 

El = .5^7 NO ~ WO + < tS~- 

A m+1 

Since in this case A m = m 4 , we have 

log Eg a i = ci - 4a i log(m + 1), 


and 

log Ei = C 2 — 4^2 log(m + 1). 

Thus, a log-log plot of the error in terms of the wave number will easily determine 
the rate of convergence. In Figure 2, we have plotted the rate of change of the graph 
of the log-log plot of the error in terms of wave number. The theory suggests that 
the rate of convergence for the Galerkin method a x is less than or almost equal to 
9/8 and for the nonlinear Galerkin method a 2 is less than or almost equal to 15/8. 
The results plotted in Figure 2 show that the Galerkin calculation asymptotes at a 
value of 1.1, whereas the nonlinear Galerkin asymptotes at a value of 1.92, which is 
in a good agreement with the theory. 

Now we turn to the case when the force is in the Gevrey class (real analytic). 
We consider the KS equation with zero forcing. With the help of the software 
package AUTO, we start the calculation with initial data on the unstable manifold 
of a periodic orbit, which again is on the attractor. The solution converges to a 
steady state as time goes to infinity. Thus, this trajectory is contained in the global 
attractor, and the theory presented in Section 3 holds for this trajectory. We first 
compute this trajectory using 100 Fourier modes, which we will consider as our 
“exact” solution. The L 2 norm of this solution vs. time is shown in Figure 3. We 
integrate this trajectory out to 1.3 time units and record the solution every .01 time 
units. 

We expect from the theory outlined in Section 3 that 


C\t 


-<*\ 


A 


1/4 
m + 1 


'm+ 1 


E 9°‘ ■= 0 max 3 M0 - MOlU* < 
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FIGURE 1. L 2 norm of the solution vs. time for the KS eq. of a Galerkin scheme 
with 100 modes forced with / = V.7 s * n (i x )- 



Log of the Number of Modes 


Figure 2. Rate of change of the log of the accuracyxlO vs. bg of the nu , caber 
f modes. The Galerkin asymptotes at 4.4 and the nonlinear Galerkin at 7.7. Th 
heory suggests that the Galerkin should asymptote near 4.5 and the nonlinear 
Salerkin near 7.5. • , Gelerkin; - , FMT. 
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C2C 


El := o“o.3 11“^) “ (*"(*) + ^n*l(*m(0))IU* < A 7 /4 


-**LXx 


m+1 


for some U\ > 0. Since A m +i = (m + l) 4 , we have 

log Egat = -Ox (m + 1) + log(ciA“+i) 


and 

log -^i = — <Tj(m + 1) + log(c 2 A“Y 1 4 ). 

Thus, a log- linear plot of the error versus the number of modes should be nearly 
linear. This is confirmed in Figure 4. In fact, the two lines are parallel. That is, 
they have the same exponential rate of convergence (same <T \ ). In addition, notice 
that the nonlinear Galerkin method still exhibits an algebraic improvement over the 
standard Galerkin. This is manifested by the fact that the graph for the nonlinear 
Galerkin is below the graph for the Galerkin. 

Similar considerations apply where the compatibility of the solutions with the 
basis functions is an issue. Let us recall how this can come about. Consider Burgers’ 
equation forced on the boundary 


du d 2 u du 
dt~ !h* + U !h = °’ 


(13) 


u(0,t) — 1 u(l,f) = 0. 

To formulate this problem in the same setting as Equation (3), one can set v = 
u — (1 — x) to obtain 


dv 

dt 


d 2 v dv dv 


— v = l — x, 


v(0, t) — 0 v(l,f) = 0. 

To apply the nonlinear Galerkin method, we must first expand the forcing term, 
here / = 1 — x, in terms of the eigenfunctions of the linear dissipative operator 
= sin(j7rx). We find 


/ = 1 — x = 



It is easily seen that, in terms of the spectral method based on the eigenfunctions 
of the Laplacian, the forcing term is only L 2 ((0, l)),(o = 0). We therefore expect 
the nonlinear Galerkin method to be significantly more accurate in this case. 
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Figure 3. L 2 norm of “exact” solution vs. time for the KS equation with zero 
forcing with initial data on an unstable manifold. 



Number of Modes 


FIGURE 4. Log of the error xlO 13 vs. the number of modes for the KS equation 
with zero forcing. Notice that the rate of convergence of both schemes is exponential. 
• , Gelerkin; ■ , FMT 
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3* Future plans 

The goal remains to understand what aspects of the long time behavior of infinite- 
dimensional PDEs axe retained by their finite-dimensional approximations. Perhaps 
the ultimate test of an approximation will remain how well it predicts nature. How- 
ever, there are a growing number of applications where one does not know a priori 
nor is there a way to test experimentally the behavior of some systems. Therefore, 
assuming the dynamics of the PDEs accurately reflects the physical phenomenon it 
is meant to depict, we hope understanding the behavior of approximating schemes 
of these PDEs will prove valuable in the future. 
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